Angular  distribution  of  Ly a  resonant  photons  emergent  from  optically  thick 

medium 

Yang  Yang1,  Ishani  Roy2,  Chi-Wang  Shu1  and  Li-Zhi  Fang3 

ABSTRACT 

We  investigate  the  angular  distribution  of  Ly  a  photons  transferring  in  or  emergent 
from  an  optically  thick  medium.  Since  the  evolutions  of  specific  intensity  I  in  the  fre¬ 
quency  space  and  the  angular  space  are  coupled  from  each  other,  we  first  develop  the 
WENO  numerical  solver  in  order  to  find  the  time-dependent  solutions  of  the  integro- 
differential  equation  of  I  in  the  space  of  frequency  and  angular  simultaneously.  We 
show  first  that  the  solutions  with  the  Eddington  approximation,  which  assume  I  to  be 
linearly  dependent  on  the  angular  variable  //,  yield  similar  frequency  profiles  of  the 
photon  flux  as  that  without  the  Eddington  approximation.  However,  the  solutions  of 
the  //  distribution  evolution  are  significantly  different  from  that  given  by  Eddington 
approximation.  First,  the  angular  distribution  of  I  are  found  to  be  substantially  de¬ 
pendent  on  the  frequency  of  photons.  For  photons  with  the  resonant  frequency  vo,  / 
contains  only  a  linear  term  of  //.  For  photons  with  frequency  at  the  double  peaks  of  the 
flux,  the  //-distribution  is  highly  anisotropic,  in  which  most  photons  are  in  the  direc¬ 
tion  of  radial  forward.  Moreover,  either  at  v0  or  at  the  double  peaks,  the  //  distributions 
actually  are  independent  of  the  initial  //  distribution  of  photons  of  the  source.  This  is 
because  the  photons  with  frequency  either  of  Vo  or  of  the  double  peaks  have  undergone 
the  process  of  forgetting  their  initial  conditions  due  to  the  resonant  scattering.  We  also 
show  that  the  optically  thick  medium  is  a  collimator  of  photons  at  the  double  peaks. 
Photons  from  the  double  peaks  form  a  forward  beam  with  very  small  spread  angle. 
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1.  Introduction 


The  transfer  of  Ly a  photons  in  an  optically  thick  medium  consisting  of  neutral  hydrogen 
atoms  is  fundamentally  important,  as  Lycr  photon  is  an  effective  tool  to  study  the  physics  of  lumi¬ 
nous  objects  at  various  epoches  of  the  universe.  It  includes  Ly  a  emitters,  Lya  blob,  damped  Ly  a 
system,  Ly  a  forest,  fluorescent  Lya  emission,  star-forming  galaxies,  quasars  at  high  redshifts  as 
well  as  optical  afterglow  of  gamma  ray  bursts  (Haiman  et  al.  2000;  Fardal  et  al.  2001;  Dijkstra  & 
Loeb  2009;  Latif  et  al.  2011).  Lya  photons  emergent  from  an  optically  thick  halo  surrounding  a 
source  of  the  Lya  photon  would  be  useful  to  constrain  the  mass  density,  velocity,  temperature  and 
the  fraction  of  neutral  hydrogen  of  the  optically  thick  halo. 

The  radiative  transfer  of  resonant  photons  has  been  extensively  studied  analytically  and  nu¬ 
merically  for  more  than  half  a  century.  An  early  effort  (Adams  et  al.  1971)  only  focused  on  the 
numerical  approximation  of  the  redistribution  function  of  resonant  scattering,  and  no  solution  of 
the  integro-differential  equation  of  the  radiative  transfer  has  been  found.  The  first  analytical  solu¬ 
tion  of  the  integro-differential  equation  is  given  by  Field  (1958)  for  the  case  of  both  medium  and 
source  to  be  uniformly  distributed  in  the  whole  space.  Analytical  solutions  of  the  frequency  profile 
of  photons  emergent  from  optically  thick  halo  are  also  found  based  on  the  Fokker-Planck  (P-F) 
approximation  of  the  integro-differential  equation  (Harrington  1973;  Neufeld  1990;  Dijkstra  et  al. 
2006).  Monte  Carlo  (MC)  simulations  are  also  popular  in  solving  the  transfer  of  resonant  photons 
(e.g.  Loeb  &  Rybicki  1999;  Zheng  &  Miralda-Escude  2002;  Tasitsiomi  2006;  Verhamme  et  al. 
2006;  Laursen  &  Sommer-Larsen  2007;  Pierleoni  et  al.  2009;  Xu  &  Wu  2010;  Xu  et  al.  2011). 

Nevertheless,  many  important  topics  cannot  be  seen  with  the  above-mentioned  solutions.  Be¬ 
sides  the  Field’s  analytical  solution,  all  others  are  time-independent,  and  therefore,  they  cannot 
even  be  used  to  describe  the  formation  and  evolution  of  the  Wouthuysen-Field  (W-F)  local  ther- 
malization  of  the  Ly  a  photon  frequency  distribution  (Wouthuysen  1952;  Field  1958,  1959).  The 
rich  features  of  the  Ly  a  photon  transfer  referring  to  the  W-F  local  thermalization  are  fully  missed. 
The  P-F  equation  is  based  on  the  Eddington  approximation,  which  assumes  that  the  radiation  in¬ 
tensity  is  a  linear  function  of  angular  (direction)  variable.  The  solutions  of  the  P-F  equation  do  not 
provide  the  information  of  the  evolution  of  the  angular  distribution  of  Ly  a  photons. 

Recently,  a  solver  of  the  radiative  transfer  of  resonant  photon  has  been  developed  based  on 
the  weighted  essentially  non-oscillatory  scheme  (WENO),  which  is  a  good  numerical  solver  of  the 
hydrodynamic  equations  and  the  kinetic  equations  (Jiang  &  Shu  1996).  With  the  WENO  solver, 
many  interesting  features  of  Lycr  resonant  photon  transfer  have  been  revealed.  It  shows  that  the 
double  peaked  frequency  profile  of  the  Lya  photon  emergent  from  an  optically  thick  medium 
generally  does  not  follow  the  time-independent  solutions  of  the  P-F  equation  (Fang  2009;  Roy 
et  al.  2009a,  2009b,  2009c,  2010).  The  solver  has  also  been  used  to  calculate  time-dependent 
features:  1)  the  timescale  of  the  formation  of  the  W-F  effect;  2)  the  light-curve  of  a  flash  source 
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surrounded  by  optically  thick  halo;  3)  the  dependence  of  the  Lya  photon  escape  coefficient  on  the 
redshifts  of  sources;  4)  the  effects  of  dust  on  the  double-peak  profile  of  the  emergent  Lyo-  photons 
(Roy  et  al.  2009a,  2009b,  2009c,  2010;  Yang  et  al.  2011;  Xu  et  al.  2011).  These  features  have  not 
been  seen  with  previous  methods. 

This  paper  will  study  the  angular  distribution  of  Lya  photon  transferring  in  an  optically  thick 
medium.  As  mentioned  above,  the  evolution  of  the  angular  distribution  is  completely  ignored  in  all 
the  methods  based  on  the  Eddington  approximation.  The  evolution  of  angular  distribution  actually 
is  significant.  In  a  thermalized  or  statistical  equilibrium  state,  the  angular  distribution  of  photons 
should  be  isotropic,  regardless  of  the  initial  angular  distribution.  Therefore,  one  can  expect  that  the 
angular  distribution  of  Lycr  photons  with  resonant  frequency  v0  should  be  isotropic.  On  the  other 
hand,  the  angular  distribution  of  photons  with  frequency  different  from  v0  might  be  anisotropic,  as 
those  photons  are  not  involved  in  the  evolution  of  thermalization  or  statistical  equilibrium.  Con¬ 
sequently,  the  angular  distributions  of  Lycr  resonant  photons  from  optically  thick  medium  should 
be  frequency-dependent.  It  definitely  cannot  be  described  by  the  Eddington  approximation.  The 
evolution  of  the  angular  distribution  of  resonant  photons  is  not  trivial.  We  still  use  the  WENO 
solver.  In  this  case,  we  should  first  develop  the  WENO  solver  to  be  able  to  simultaneously  solve 
the  photon  transfer  in  both  frequency-  and  angular  spaces. 

The  paper  is  organized  as  follows.  Section  2  presents  the  basic  features  of  Lya  photon  transfer 
given  by  the  WENO  solver.  In  section  3,  the  precision  of  the  Eddington  approximation  will  be  stud¬ 
ied  via  a  comparison  between  solutions  with  and  without  the  Eddington  approximation.  Section 
4  presents  the  results  of  the  evolution  of  the  angular  distribution,  especially,  on  the  frequency- 
dependence,  source-dependence  and  position-dependence  of  the  Lycr  photon  angular  distribution. 
The  discussion  and  conclusion  are  given  in  Section  5.  Mathematical  details  of  the  WENO  algo¬ 
rithm  on  the  radiative  transfer  equation  are  given  in  the  Appendix. 


2.  WENO  solver  of  transfer  equations  of  resonant  photons 

The  WENO  solver,  some  details  of  which  being  given  in  the  Appendix,  is  to  solve  the  follow¬ 
ing  radiative  transfer  equation  of  Ly a  resonant  photon  in  a  spherical  symmetric  medium  containing 
neutral  HI. 

dl  dl  (1  -n2)dl  _  dl  _ 
dr)  ^  dr  r  d/u  ^  dx 

-(p(x;a)I  +  J'  K(x,  x'\ a)I(rj,  r,  x' ,n')dx'dn' /2  +  S,  (1) 

where  I(t,rp,x,/u )  is  the  specific  intensity,  which  is  a  function  of  time  /,  radial  coordinate  rp, 
frequency  x  and  /./  =  cos  0,  the  direction  angle  with  respect  to  the  radial  vector  r.  S  (//,  rp,  x,/r)  is 
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the  source  of  photons. 

In  eq.(l),  we  use  the  dimensionless  time  ?/  defined  as  tj  =  cnmcr0t  and  the  dimensionless  radial 
coordinate  r  defined  as  r  =  nmcr{)rp,  where  nm  is  the  number  density  of  HI,  and  cr0/^1/2  is  the  cross 
section  of  HI  resonant  scattering  of  Lycr  photons  at  resonant  frequency  vo  =  2.46  x  1015  s-1.  That 
is,  T]  and  r  are,  respectively,  in  the  units  of  mean  free  flight-time  and  mean  free  path  of  photon  v0 
with  respect  to  the  resonant  scattering  without  dust  scattering  and  absorption.  Without  resonant 
scattering,  a  signal  propagates  in  the  radial  direction  with  the  speed  of  light,  the  orbit  of  the  signal 
is  then  r  =  tj  +  const. 


(p{x,  a)  is  the  normalized  Voigt  profile  (Hummer  1965)  given  as 


1 p(x ,  a) 


a 

n 3/2 


dy 


{x  -y)2  +  a 


2 • 


(2) 


As  usual,  the  photon  frequency  v  in  eq.(2)  is  described  by  the  dimensionless  frequency  x  =  (v  - 
vo)/Avd,  and  Avr)  =  vo(vT/c )  =  1.06  x  10'  '(T / 104)1/2  Hz  is  the  Doppler  broadening  by  the  thermal 
motion  vT  =  yj2kBT/mm,  T  being  the  gas  temperature  of  the  halo.  The  parameter  a  in  eq.(2)  is  the 
ratio  of  the  natural  to  the  Doppler  broadening.  For  the  Lyo-  line,  a  =  4.7  x  HT4(7 / 1 04)~ 1  /2.  The 
optical  depth  of  Lya  photons  with  respect  to  HI  resonant  scattering  is  drs(x)  =  nmcr(x)drp,  where 
cr(.r)  =  o-q4>(x,  a)  is  the  cross  section  of  scattering  at  v,  and  therefore,  the  dimensionless  size  of  the 
halo  R  is  equal  to  the  optical  depth  r0  =  nm(r{)R. 


The  re-distribution  function  R(x,x';a)  of  eq.(2)  gives  the  probability  of  a  photon  absorbed 
at  the  frequency  x',  and  re-emitted  at  the  frequency  x.  It  depends  on  the  details  of  the  scattering 
(Henyey  &  Greestein  1941;  Hummer  1962;  Hummer  1969).  If  we  consider  coherent  scattering 
without  recoil,  the  re-distribution  function  with  the  Voigt  profile  is 


R(x,x';a)  = 

1  r°° 

^  X-v 


\x-x'\/2 


tan 


+  u 


-  tan 


du. 


(3) 


where  xmin  =  minG,  x')  and  ,rmax  =  maxCr,  x').  In  the  case  of  a  =  0,  i.e.  considering  only  the 
Doppler  broadening,  the  re-distribution  function  is 


«(„')  =  -erfc[max(W,M)]. 


(4) 


The  re-distribution  function  of  equation  (4)  is  normalized  as  R(x,  x')dx'  =  (Mx,  0)  =  n~x,2e~x\ 
With  this  normalization,  the  total  number  of  photons  is  conserved  in  the  evolution  described  by 
equation  (3).  That  is,  the  destruction  processes  of  Lya  photons,  such  as  the  two-photon  process 
(Spitzer  &  Greenstein  1951;  Osterbrock  1962),  are  ignored  in  equation  (3).  The  recoil  of  atoms  is 
not  considered  in  equation  (3)  or  (4).  The  dust  absorption  is  also  ignored. 
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For  Dipole  scattering,  eq.(4)  should  be  replaced  by  (Hummer  1962) 

3(1 

R{x,  x')  =  -  <  -erfc[max(|jc|,  |jc'|)][3  +  2{x2  +  x'2)  +  4x2x'2]  (5) 

-  — — — p— —  exp[-max(|.r|,  |.r'|)][2|min(|x|,  |.r'|)|2  +  1] 

X/7T 

which  is  still  /./-independent. 

In  eq.  (1),  the  term  with  the  parameter  y  is  due  to  the  expansion  of  the  universe.  If  «hi  is 
equal  to  the  mean  of  the  number  density  of  cosmic  hydrogen,  we  have  y  =  tq1p,  and  tqp  is  the 
Gunn-Peterson  optical  depth.  Since  the  Gunn-Peterson  optical  depth  is  of  the  order  of  106  at  high 
redshift  (e.g.  Roy  et  al.  2009c),  the  parameter  y  is  of  the  order  of  10  5  -  1 0-6 .  Therefore,  if  the 
optical  depth  of  halos  is  equal  to  or  less  than  106,  the  term  with  y  of  eq.(l)  can  be  ignored. 


2.1.  Test  with  Field’s  analytical  solution 


The  WENO  solver  used  in  this  paper  is  different  from  the  previous  version  (Roy  et  al.  2009a, 
2009b,  2009c).  In  the  previous  version,  the  evolution  of  /  in  the  //- space  is  eliminated  by  the 
Eddington  approximation,  while  the  current  solver  treats  the  evolution  of  I  in  the  //-space  governed 
by  eq.(l). 


We  first  test  the  WENO  solver  with  analytical  solutions.  Assuming  that  the  specific  intensity 
and  source  S  are  homogeneous  in  the  r  and  //  space,  i.e.  /(//,  r,  x,jj)  is  independent  of  variables  r 
and  //.  Eq.(l)  becomes 


where 


dJ  dJ  r  , 

— - y—  =  —<p(x)I  +  R(x,  x  x  )dx  +  S, 

OT]  OX  J 


J(r h  x) 


u 


I(r/,r,x,n)dn. 


(6) 

(7) 


Take  y  =  0,  Doppler  parameter  a  =  0,  the  source  S  =  <p( x)  =  n  1/2 e  x~ ,  and  the  initial  radiative  field 
I(x,  rj  -  0)  =  0.  The  time-dependent  solution  of  eq.(6)  is  (Field  1958,  Rybicki  &  DelTAntonio 
1994) 


J(x,tj )  =  n  1/2[1  -  exp (-Tje  *")] 


+ 


ew~[l  -  (1  +  rje  "")exp {-rje  w~)]erf (w)dw. 


(8) 


Our  solver  is  to  directly  find  the  solution  I  from  eq.(l).  One  can  then  give  J  via  eq.(7).  It  is 
interesting  to  see  whether  the  solution  eq.(8)  can  be  reproduced,  if  we  also  assume  that  the  source 
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S  in  eq.(l)  is  spatially  homogeneous,  but  /^-dependent,  i.e.  S  =  Q(ju)(p(x)  =  Q(/i)n~'/2e~x~  where 
0(//)  describes  the  angular  distribution  of  photons  from  the  source.  We  consider  the  source  as 
follows, 

„  J  2(n  +  l)n"n~l,2e~xl ,  0  <  fi  <  l, 

{  0,  -l<yU<0, 

where  n  is  taken  to  be  a  positive  integer.  Obviously,  the  source  is  isotropic  when  n  =  0,  while  it  is 
anisotropic  if  n  =£  0.  The  larger  the  n,  the  stronger  the  emission  in  the  direction  n  =  1.  The  factor 
2(n  +  1)  is  for  normalization:  2  f  2 (n  +  l)iund/u  =  1. 

The  numerical  results  with  sources  n  =  0,  4  and  6  are  shown  in  Fig.  1.  It  is  expected  that  the 


Fig.  1. —  The  WENO  numerical  solutions  (solid  lines)  of  eq.(l)  assuming  the  sources  S  is  a) 
S  =  n~x,1e~x~  for  all  n  (left);  b)  S  =  I0fi4n~1/2e~x'  (middle);  and  c)  S  =  14 /j.6n~l/2e^x~  (right)  for 
/u  >  0  and  S=0  for  /i  <  0.  The  Field’s  analytical  solution  is  shown  with  dot  lines. 

numerical  solution  of  n  =  0  (the  left  panel  of  Fig.l)  should  follow  the  analytical  solution  eq.(8) 
well,  as  the  isotropic  source  is  the  same  as  that  used  to  find  the  analytical  solution. 

It  is  interesting  to  see  that  the  WENO  solutions  of  n  =  4  (middle  panel)  and  n  =  6  (right 
panel)  also  follow  the  analytical  solution  eq.(8)  well.  It  seems  to  indicate  that  the  evolution  of  the 
frequency  space  is  independent  of  that  of  the  yu- space. 


2.2.  Time  scale  of  the  statistical  equilibrium  of  the  angular  distribution 

A  remarkable  feature  of  the  solutions  of  Fig.  1  is  to  show  a  flat  plateau  in  the  range  |jc|  <  2  at 
time  77  >  100.  The  flat  plateau  is  caused  by  the  Wouthuysen-Field  local  thermalization  of  frequency 
distribution  of  resonant  photon  (Wouthuysen  1952;  Field  1958,  1959).  The  flat  plateau  actually  is 
the  Boltzmann  statistical  equilibrium  distribution  around  x  =  0  when  the  atomic  mass  is  infinite. 
If  the  mass  is  finite,  i.e.  considering  the  recoil  in  the  re-distribution  functions  (3)  or  (4),  the  flat 
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plateau  will  become  e~2bx,  where  b  =  hv{)/mvTc.  This  is  the  local  Boltzmann  distribution  required 
by  the  Wouthuysen-Field  effect  (Roy  et  al.  2009b).  The  resonant  scattering  between  photons  and 
HI  atoms  leads  to  the  Boltzmann  distribution  of  the  photon  frequency  distribution  around  x  =  0 
with  the  temperature  equal  to  that  of  HI  atoms. 
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Fig.  2. —  The  WENO  numerical  solutions  of  angular  distributions  from  eq.(l)  at  x  =  0,  assuming 
the  sources  S  are  a)  S  =  n~1/2e~x2  for  all  //  (left)  and  b)  S  =  10/j4n~l/2e~x2  (middle);  and  c) 
S  =  l4/j6n~1/2e~x2  (right)  for  //  >  0  and  S=0  for  //  <  0. 


When  resonant  photons  undergo  the  local  thermalization  in  the  frequency  space,  the  angular 
distribution  should  undergo  the  evolution  of  approaching  statistical  equilibrium.  The  anisotropic  //- 
distributions  have  to  evolve  to  isotropic  (statistical  equilibrium).  We  calculate  all  the  //-distributions 
at  the  time  rj  corresponding  to  the  three  panels  of  Fig.  1.  The  result  is  plotted  in  Fig.  2.  The  //- 
distribution  of  left  panel  is  always  isotropic.  This  is  expected  as  the  source  is  isotropic,  which  is 
already  in  the  state  of  statistical  equilibrium. 

The  middle  and  right  panels  of  Fig.  2  show  the  evolution  of  an  anisotropic  /./-distribution 
eq.(9)  to  isotropic.  The  time  scale  of  approaching  isotropic  distribution  seems  to  be  independent 
of  the  anisotropy  of  sources.  It  is  always  equal  to  about  77  ~  100  for  both  n  =  4  and  n  =  6,  i.e.  the 
//-distribution  will  become  isotropic  after  100  times  of  resonant  scattering.  This  time  scale  is  about 
the  same  as  that  of  the  W-F  thermalization  (Fig.  1).  Therefore,  the  thermalization  in  the  frequency 
space  and  the  isotropic  distribution  in  the  //-space  are  realized  at  about  the  same  time. 
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3.  Precision  of  the  Eddington  approximation 
3.1.  Equations  of  the  Eddington  approximation 


We  now  consider  the  transfer  of  Ly a  photons  in  a  spherical  halo  with  an  optical  source  at  its 
center.  The  halo  is  assumed  to  consist  of  uniformly  distributed  HI  gas  with  number  density  nH i- 
The  /./-dependence  of  the  specific  intensity  I  can  generally  be  expressed  by  a  Legendre  expansion 
/(//,  r,  x,  jj)  =  X/  h^h  r,  x)Pi(ju).  The  Eddington  approximation  is  to  take  only  the  first  two  terms  of 
the  Legendre  expansion,  and  drop  all  terms  of  l  >  2  (e.g.  Rybicki  &  Lightman,  1979).  That  is 

I(ji,  r,  x,  //)  -  J{ tj,  r,  x)  +  3fxF(rj,  r,  x),  (10) 


where 

1  f+1  1  r+1 

J(i h  r,x)  =  -  J  I(ji,  r,  x,  / u)dfi ,  F(tj,  r,x)=  -  J  /uHi),  r,  x,  iu)d/u. 

They  are,  respectively,  the  angularly  averaged  specific  intensity  and  flux. 


(11) 


Defining  j  =  r2J  and  /  =  r2F ,  eq.(l)  yields  the  equations  of  j  and  /  as 


di  +  df 

drj  dr 


-<P(x\  a)j  + 


r 

*R(x,  x\  a)  jdx  +  y 


dj  2 

T-  +  r  S, 
ox 


df_  +  1^7  _  2 1 
drj  3  dr  3  r 


df 

~<f>(x;  a)f  +  y— . 

dx 


(12) 

(13) 


The  mean  specific  intensity  j{rj,  r,  x )  describes  the  x  photons  trapped  in  the  position  r  at  time  // 
by  the  resonant  scattering,  while  the  flux  f(rj,  r,  x)  describes  the  photons  in  transit.  One  can  test 
the  precision  of  the  Eddington  approximation  by  a  comparison  of  the  solution  of  eq.(l)  without 
Legendre  expansion  with  that  of  the  Eddington  approximation. 


3.2.  Profiles  of  j  and  / 

For  spherical  halo  with  a  central  source,  the  term  S  of  eq.(l)  can  be  replaced  by  a  boundary 
condition  of  I(jj,r,x,iu )  at  r  =  0.  If  the  angular  distribution  of  photon  is  independent  of  photon’s 
frequency,  we  have  generally 


r2I(rj,r,x,n)\r-,o  =  S0T(T])®(jF)(p(x).  (14) 

where  the  functions  T (//),  ©(//),  and  0(a)  describe,  respectively,  the  time-dependence,  angular-  and 
frequency-distributions  of  photons  of  the  source.  The  factor  S  0  gives  the  intensity  of  the  source. 
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In  this  case,  the  source  of  eq.(12)  can  be  replaced  by  a  boundary  condition  at  r  =  0  as 

7(77,0,*)  =  0,  f(r],0,x)  =  S 0T(ri)(p(x)^  J' pQ(p)dp.  (15) 

On  the  outside  of  the  halo,  r  >  R,  no  photons  propagate  in  the  direction  p  <  0.  The  boundary 
condition  at  r  -  R  of  eq.(l)  should  be 

I(rj,R,x,p)  =  0,  p  <  0.  (16) 


/■» —  1 

For  eq.(12),  we  have  J  pI(jj,R,  x,p)dp  =  0  (Unno  1955),  the  boundary  condition  is  then 


j(jl,R,x)  =  2  f(rj,R,x). 


If  the  source  starts  to  emit  photon  at  t  =  0,  the  initial  condition  should  be 

l(0,r,x,p )  =  0, 


for  eq.(l),  and 

7(0,  r,  x)  =  /( 0,  r,  x)  =  0, 

for  eq.(12). 

We  solve  eq.(l)  by  taking  the  boundary  condition  at  r  =  0  to  be 


r2I(ri,r,x,p)\r^o  = 


p>0, 

0,  p  <  0. 


(18) 

(19) 

(20) 


With  eq.(20)  one  can  find  I  from  eq.(l),  and  then  find  j  and  /  via  eqs.(l  1).  The  results  are  given 
in  the  right  panels  of  Fig.  3,  which  shows  the  well-known  double  peaks  at  |jc|  —  2. 


With  eq.(20),  the  corresponded  boundary  condition  of  eq.(15)  is 

j(R,  r  =  0,  x)  =  0,  f(j],  r  =  0,  x)  =  n~'/2e~x\  (21) 


The  profiles  of  j  and  /  given  by  the  solutions  of  eqs.(l  1)  and  (12)  and  condition  eq.(21)  are  shown 
in  the  left  panels  of  Fig. 3.  The  profiles  on  the  left  and  right  panels  in  Fig. 3  are  about  the  same, 
indicating  that  the  Eddington  approximation  is  a  good  one  for  this  case. 

To  replace  the  factor  6 p  of  eq.(20)  with  Q(p)  =  1 6//’,  we  re-do  the  solution  of  j  and  /  with 
eq.(l).  The  results  are  given  in  Fig. 4,  which  still  shows  the  same  shape  of  the  profiles  as  that  of 
the  left  panels  of  Fig. 3.  That  is,  the  profiles  of  j  and  /  are  not  affected  by  the  angular  distribution 
of  photons  from  the  source.  It  is  probably  because  the  /./-distribution  quickly  evolves  into  the 
statistical  equilibrium  state,  the  initial  anisotropy  of  the  p  distribution  is  forgotten. 
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Fig.  3. —  A  comparison  of  the  solutions  j  (top)  and  /  (bottom)  with  the  Eddington  approxima¬ 
tion  (left)  and  the  solutions  of  j  and  /  without  the  Eddington  approximation  (right).  Relevant 
parameters  are  R  =  102,  and  a  =  10“3. 

Nevertheless,  Figs.  3  and  4  do  show  small  differences  between  the  solutions  with  and  without 
the  Eddington  approximation,  even  though  both  solutions  are  given  by  the  same  source.  The 
difference  comes  from  the  contribution  of  the  terms  of  /  >  2  in  the  Legendre  expansion.  The 
difference  between  the  profiles  with  and  without  the  Eddington  approximation  becomes  smaller 
when  the  time  r/  is  larger.  It  is  because  larger  //  corresponds  to  larger  optical  depth.  The  Eddington 
approximation  generally  is  good  for  optically  thick  medium. 
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Fig.  4. —  The  same  as  the  right  panels  of  Fig.  3,  but  the  //-distribution  6 //  of  the  source  eq.(20)  is 
replaced  by  16/A 


4.  Angular  distributions 
4.1.  Frequency  dependence 

Although  the  Eddington  approximation  is  acceptable  to  calculate  the  profile  of  Ly a  photons  in 
the  frequency  space,  it  would  fail  in  the  /./-space.  The  result  of  Fig.  2  shows  that  the  //-distribution 
is  isotropic  at  frequency  v0.  On  the  other  hand,  the  //-distribution  will  no  longer  be  isotropic  at 
frequency  |.v|  >  2,  because  photons  of  |.i'|  >  2  have  not  undergone  enough  number  of  scattering. 
Consequently,  the  angular  distribution  of  photons  emergent  from  optically  thick  halo  should  be 
frequency(energy)-dependent. 

We  calculate  the  //-distribution  of  photons  from  halo  with  R  =  500  with  the  central  source 
given  by  eq.(20),  i.e.  photons  from  the  source  can  be  described  by  the  Eddington  approximation 
eq.(10).  The  result  is  shown  in  Fig. 5.  The  //  distributions  at  frequencies  x  =  0  (top  left  panel  of 
Fig.  5)  and  0.8  (top  right)  basically  are  straight  lines  in  the  whole  range  -1  <  //  <  1.  That  is,  I  can 
be  described  by  the  Eddington  approximation  eq.(10). 

At  x  =  1.6  (bottom  left  panel  of  Fig.  5),  the  //-distribution  starts  to  deviate  from  a  straight 
line,  i.e.  deviating  from  an  Eddington  approximation.  At  x  =  2.4  (bottom  right),  the  //-distribution 
shows  a  very  sharp  spike  at  //  =  1.  That  is,  the  angular  distribution  of  photons  with  frequency 
at  the  two  peaks  (Fig.  3)  is  significantly  different  from  isotropic,  but  is  dominated  by  photons  of 
//  =  1.  This  result  is  consistent  with  the  “single  shot  picture”  (Adams  1975;  Bonilha  et  al.  1979), 
in  which  photons  with  frequency  |jc|  <  2  mainly  undergo  a  diffusion  in  the  frequency  space;  once 
a  photon  diffuses  to  |jc|  >  2,  it  will  take  “single  longest  excursion”  to  leave  for  outside  of  the  halo. 


-12- 


Fig.  5. —  The  fj  distribution  of  photons  emergent  from  a  halo  with  radius  R  =  500.  The  frequencies 
are  x  =  0.0  (top  left);  0.8  (top  right);  1.6  (bottom  left)  and  2.4  (bottom  right).  The  relevant 
parameters  of  the  calculation  are  r)  =  1.2  x  104,  y  =  0,  and  a  =  10~3. 

Therefore,  the  two  peaks  of  the  flux  /  at  frequency  x+  -  ±(2  -  3)  are  dominated  by  photons  from 
“single  longest  excursion”  photons,  of  which  n  ~  1. 


4.2.  Dependence  of  the  initial  anisotropy 


The  source  of  Fig.  5  given  by  eq.(20)  has  ®{/u)  =  6/u  (ju  >  0),  which  is  linear  of  /r.  We  now 
consider  sources  with  higher  anisotropy  with  ®(ju)  given  by 


®(M) 


2 (n  +  2)/r",  0  <  n  <  1, 

0,  /u  <  0. 


(22) 


When  the  integer  n  is  large,  0(/r)  is  similar  to  a  6  function  6(ju  -  1),  i.e.  most  photons  are  in  the 
direction  //  =  1 . 
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Fig.  6. —  The  //-distribution  of  halo  with  radius  R  =100  and  source  eq.(22)  with  n- 1,  2,  4,  and 
6.  The  frequencies  are  taken  to  be  x  =  0  (top  left),  0.8  (top  right),  1.6  (bottom  left),  2.4  (bottom 
right).  The  parameters  of  the  halos  are  tj  =  3500,  y  =  0,  and  a  =  10”3. 

We  repeat  the  calculation  of  Fig. 5,  but  using  the  source  eq.(22)  with  n  =  1,  2,  4  and  6.  The 
result  is  plotted  in  Fig.  6.  It  is  interesting  to  see  that  the  //-distributions  are  independent  of  n,  but 
only  depend  on  x.  It  is  easy  to  explain  the  //-independence  of  the  two  top  panels  of  Fig. 6,  both 
of  which  have  frequency  |jc|  <  2.  In  this  frequency  range,  the  evolution  of  the  specific  intensity 
I  is  governed  by  the  local  thermalization  of  v-spacc  and  entropy  increasing  of  //-space.  These 
processes  lead  to  the  Boltzmann  distribution  in  the  energy  space  and  isotropic  distribution  in  the 
angular  space,  regardless  of  the  initial  distributions  in  either  frequency-  or  angular  spaces.  In 
other  words,  the  initial  distribution  is  forgotten  during  the  local  thermalization  and  approaching 
statistical  equilibrium. 

However,  the  mechanism  of  the  local  thermalization  and  approaching  isotropic  distribution 
seems  to  be  unable  to  explain  why  the  two  bottom  panels  of  Fig.  6.  also  show  //-independence. 
The  //-distribution  of  the  two  bottom  panels  of  Fig.  6  are  highly  anisotropic.  Therefore,  they  do  not 
have  to  be  the  result  of  the  local  thermalization  and  approaching  statistical  equilibrium.  Why  do 
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they  also  show  the  behavior  of  forgetting  the  initial  angular  distributions?  The  reason  is  as  follows. 
In  the  first  phase  of  resonant  photon  evolution,  Lycr  photons  are  trapped  in  the  range  of  |jt|  <  2 
within  the  time  scales  of  a  few  tens  or  hundred  scattering  (Roy  et  al.  2009c).  The  trapped  photons 
have  already  forgotten  their  initial  state.  On  the  other  hand,  photons  with  |jc|  >  2  mostly  come 
from  the  diffusion  of  trapped  photons  from  |jc|  <  2  to  |jc|  >  2  (Roy  et  al.  2009b).  Thus,  all  photons 
of  |.v|  >  2  emergent  from  the  optically  thick  halo  essentially  have  the  same  initial  condition,  given 
by  the  |jc|  space  diffusion  of  trapped  photons.  Therefore,  the  initial  distributions  before  they  are 
trapped  have  been  forgotten.  This  property  can  also  be  seen  in  Fig.  2,  in  which,  although  the 
sources  of  the  middle  and  right  panels  are  different  from  each  other,  the  behaviors  of  the  time- 
evolution  of  the  //-distribution  are  about  the  same.  This  result  also  implies  that  it  is  impossible  to 
find  the  information  of  the  distribution  of  photons  emitted  by  the  central  source. 


4.3.  Collimation  of  photons  of  the  double  peaks 

A  common  feature  of  the  two  bottom  panels  of  Fig.  6  is  to  show  a  very  sharp  spike  at  //  ~  1. 
The  photon  frequencies  of  the  two  panels  are  at  \x\  =  1.6  and  \x\  =  2.4,  corresponding  to  the 
double  peaks  of  Figs  3  and  4.  Therefore,  the  spiky  distribution  of  //  indicates  that  the  photons  with 
frequency  at  the  double  peaks  have  formed  a  forward  beam. 

In  order  to  measure  the  angular  size  of  the  //  =  1  spikes,  we  fit  the  //-distributions  of  the  two 
bottom  panels  of  Fig.  6  with  polynomials  of  //.  We  find  that  both  the  two  bottom  panels  of  Fig.  6 
can  be  well  fitted  with  polynomials  of  //  having  leading  terms  A//16  +  B/u15  +  ...,  A,  B  being  fitting 
coefficients.  The  terms  of  either  //16  or//15  are  much  sharper  than  the  central  source  eq.(22)  //''  with 
n  <  6.  Therefore,  the  radiative  transfer  at  the  double  peaks  of  frequency  space  plays  the  role  of 
forward  collimator.  It  made  the  photons  to  form  forward  beams. 

If  we  define  the  spread  angle  / 3  of  the  forward  beam  as  the  angle  of  half  intensity,  this  number 
can  be  estimated  by  sin/3  ~  (1/2)16,  and  therefore,  (3  ~  1(T4  angular  degree.  This  result  is  again 
consistent  with  the  “single  shot  picture”.  The  double  peaks  mainly  consist  of  photons  from  a  single 
shot,  which  moves  in  the  forward  direction. 


4.4.  Large  halo 

We  calculate  the  //-distribution  of  Lya  photons  in  a  halo  with  large  radius  R  =  1000,  and  the 
central  source  is  given  by  eq.(22)  and  n  =  6.  The  results  are  given  in  Fig.  7,  which  shows  the 
dependence  of  the  //  distribution  on  the  radial  variable  r  in  the  halo. 
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Fig.  7. —  The  //-distributions  at  radial  positions  r  -  100,  300,  500  of  a  halo  with  radius  R  =  1000. 
The  source  is  given  by  eq.(22)  and  n  =  6.  The  frequencies  are  taken  to  be  x  =  0  (top  left),  0.8  (top 
right),  1.6  (bottom  left),  2.4  (bottom  right).  Other  relevant  parameters  are  r\  -  1.2  x  104,  y  =  0  and 
a  =  10“3. 

Although  the  photons  from  the  source  of  n  =  6  are  highly  anisotropic,  all  the  //-distributions 
of  x  =  0.0  and  0.8  at  r  =  100,  300,  500  are  straight  lines.  That  is,  the  specific  intensity  I  can  be 
well  approximated  by  the  Eddington  approximation  eq.(10).  This  result  is  consistent  with  §4.2. 
The  slope  of  the  straight  line  is  decreasing  with  the  radial  position.  It  is  simply  because  for  a 
stable  solution,  the  photon  number  is  conserved,  and  the  flux  should  be  smaller  at  larger  r.  The 
r-dependence  of  the  //-distribution  of  |jt|  =  2.4  photons  is  also  consistent  with  the  result  of  section 
4.3:  the  larger  the  r,  the  sharper  the  //-distribution.  The  r  transfer  leads  to  the  collimation. 

The  behavior  of  r-dependence  of  the  //-distribution  at  x  =  1.6  (bottom  left  of  Fig.  7)  is  very 
different  from  that  of  x  =  0.0,  0.8,  and  2.4.  The  //  distribution  at  r  =  100  is  about  the  same  as  the 
bottom  right  panel  of  Fig.  6,  i.e.  it  has  undergone  an  evolution  of  forward  collimation,  having  a 
sharp  spike  at  //  =  1 .  However,  the  //  distribution  will  no  longer  show  a  spike  at  r  =  300  and  500. 
That  is,  the  r-dependence  of  the  bottom  left  panel  shows  two  phases.  When  r  is  equal  to  or  less 
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than  about  200,  the  r-dependence  of  the  /./  distribution  is  similar  to  the  |jc|  =  2.4  photons.  While 
when  r  >  200,  the  r-dependence  of  fj  distribution  is  similar  to  the  |jc|  =  0  and  0.8  photons.  This  is 
because  the  optical  depth  at  \x\  -  1.6  is  larger  than  \x\  =  2.4,  the  “single  shot  picture”  is  working 
well  at  r  ~  100  for  both  |jc|  =  1.6  and  2.4.  However,  at  r  >  200,  the  single  shot  picture  can  still 
work  well  for  |.r|  =  2.4,  but  not  so  well  for  |jc|  =  1.6. 


4.5.  Effect  of  anisotropic  scattering 

All  calculations  in  the  previous  sections  are  based  on  the  re-distribution  function  eq.(4),  which 
considered  only  isotropic  scattering.  If  we  consider  dipole  scattering,  it  seems  to  introduce  a  new 
factor  leading  to  anisotropic.  However,  the  re-distribution  function  of  dipole  scattering  eq.(5) 
actually  is  independent  of  /i.  It  does  not  change  any  of  the  results  mentioned  above.  The  dipole 
scattering  of  individual  atoms  may  yield  new  anisotropic  behavior.  However,  HI  atoms  are  in 
thermal  equilibrium,  and  their  distribution  is  isotropic.  The  dipole  scattering,  in  average,  does  not 
contain  any  parameter  of  specific  direction.  It  will  not  add  any  anisotropic  behavior.  Therefore,  all 
conclusions  in  the  previous  sections  should  still  hold. 


5.  Conclusion 

The  transfer  of  Ly a  resonant  photons  from  a  central  source  in  a  halo  consisting  of  HI  gener¬ 
ally  is  considered  as  a  problem  of  radiative  transfer  in  an  optically  thick  medium.  However,  the 
“optically  thick  medium”  is  true  only  when  the  frequency  of  Lycr  photons  lies  in  a  narrow  range 
|jc|  <  2.  The  cross  section  of  resonant  scattering  is  very  sensitive  to  the  photon  frequency.  It 
quickly  becomes  small  when  the  frequency  of  Lycr  photons  has  only  a  small  deviation  from  the 
range  |jc|  <  2.  For  those  photons,  the  halo  is  optically  moderate  thick,  or  even  thin.  Therefore,  in 
order  to  understand  the  transfer  of  Lycr  photons  with  frequency  around  the  resonant  peak,  we  need 
to  find  the  solutions  of  the  integro-differential  equation  (1)  available  simultaneously  in  optically 
thick  as  well  as  moderate  thick  and  even  thin  medium.  That  is,  although  halo  is  optically  thick  for 
resonant  photons,  one  should  not  treat  eq.(l)  by  using  the  condition  of  optical  thick. 

To  find  solution  of  eq.(l)  having  desired  precision  in  frequency  ranges  of  optically  thick  as 
well  as  moderately  thick,  the  algorithm  should  be  able  to  handle  the  extremely  flat  distribution 
(|jc|  <  2)  and  its  sharp  boundary  (|jc|  ~  2)  of  /.  These  features  can  be  properly  captured  by  the 
state-of-the-art  numerical  method,  WENO  scheme,  as  it  has  high  order  of  accuracy  and  good  con¬ 
vergence  in  capturing  discontinuities  as  well  as  to  be  significantly  superior  over  piecewise  smooth 
solutions  containing  discontinuities.  The  WENO  solver  has  been  shown  to  be  powerful  to  solve 
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the  integro-difFerential  equation  of  radiative  transfer  of  resonant  photons  Lycr.  In  this  paper,  we 
develop  the  WENO  algorithm  to  be  able  to  solve  the  integral-differential  equation  (1)  in  frequency 
and  angular  space  simultaneously. 

We  first  show  that  the  Eddington  approximation  can  yield  reasonable  results  of  the  frequency 
profile  of  photons  emergent  from  optically  thick  halos.  Since  the  Eddington  approximation  is  to 
assume  that  I  linearly  depends  on  /u,  all  the  physics  of  the  angular  distribution  of  Lya  photons  are 
missing.  A  cost  of  the  Eokker- Planck  equations  is  also  to  ignore  all  the  effects  of  the  evolutions  of 
angular  distribution. 

The  physics  of  the  evolution  of  angular  distribution  is  rich.  As  has  been  known,  resonant  scat¬ 
tering  makes  the  transfers  of  resonant  photons  in  the  physical  space  and  the  frequency  space  to  be 
coupled  between  each  other.  We  show,  in  this  paper,  that  the  resonant  scattering  leads  to  the  cou¬ 
pling  between  the  evolutions  of  resonant  photons  in  the  frequency  space  and  the  angular  space  as 
well.  The  evolution  of  resonant  photon  distribution  in  the  n  space  is  significantly  dependent  on  the 
frequency.  Photons  with  frequency  \x\  <  2  undergo  the  procedure  of  approaching  statistical  equi¬ 
librium,  and  their  angular  distribution  is  isotropic  after  a  few  tens  or  hundred  scattering,  regardless 
of  whether  the  initial  angular  distribution  is  isotropic.  On  the  other  hand,  the  angular  distribution  of 
photons  with  |jc|  >  2  is  substantially  anisotropic,  even  if  the  initial  angular  distribution  is  isotropic. 

An  interesting  feature  is  that  the  anisotropic  angular  distributions  at  frequency  |*|  ~  2  are 
independent  of  the  initial  angular  distributions.  Different  initial  angular  distributions  yield  the 
same  anisotropic  angular  distributions  after  a  few  tens  or  hundred  scattering.  This  is  because 
photons  at  frequency  |.v|  ~  2  are  not  directly  from  the  source,  but  come  from  the  trapped  photons 
within  |;c|  <  2,  for  which  the  initial  distributions  have  been  forgotten.  Therefore,  it  seems  to  be 
impossible  to  find  the  property  of  the  source  with  the  observed  p-distribution  of  Lya  photons  either 
in  the  range  of  |*|  <  2  or  in  |*|  >  2. 

Another  interesting  feature  of  an  optically  thick  halo  is  the  collimation  of  photons  with  fre¬ 
quency  of  the  double  peaks.  This  is  also  because  photons  trapped  in  |*|  <  2  are  thermal.  When 
the  trapped  photons  diffuse  to  |*|  >  2,  they  have  two  possible  fates.  One  is  to  get  out  of  the  holes 
by  a  single  shot  if  photons  move  forward.  If  a  photon  has  not  taken  a  single  shot,  the  resonant 
scattering  will  lead  it  to  get  back  to  the  region  of  |*|  <  2.  Therefore,  photon  transfer  in  optically 
thick  medium  is  a  collimator.  Although  photons  stored  in  an  optically  thick  halo  are  thermal,  the 
/i-distribution  is  isotropic,  the  double  peak  is  only  to  pick  up  photons  of  single  shot,  i.e.  moving 
forward. 


Acknowledgment:  This  research  is  partially  supported  by  ARO  grants  W911NF-08-1-0520  and 
W911NF-1 1-1-0091. 
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A.  Numerical  algorithm 

To  solve  equation  (1),  our  computational  domain  is  (r,x,/j)  e  [0,  rmax]  x  [jqeft,  bright]  x 
[-1,1],  where  /'max,  -qcft  and  are  chosen  such  that  the  solution  vanishes  to  zero  outside 

the  boundaries.  We  choose  mesh  sizes  with  grid  refinement  tests  to  ensure  proper  numerical  res¬ 
olution.  In  the  following,  we  describe  numerical  techniques  involved  in  our  algorithm,  including 
approximations  to  spatial  derivatives,  numerical  integration,  numerical  boundary  condition  and 
time  evolution. 


A.l.  Conservation  law 

To  perform  the  WENO  algorithm,  we  first  need  to  rewrite  the  equation  into  the  form  of  a 
conservation  law.  Noticing  the  boundary  condition  (14),  we  define  /'  =  r2/,  then  equation  (1) 
becomes 

ar  or  \d(\-n2)r  dr 

— - 1-  u  — - 1 - y -  = 

drj  dr  r 


dn  dx 

-cp(x\  d)l'  +  J'  K(x,  x';a)I'(r],  r,  x' ,n')dx'dn' /2  +  r2S . 

For  simplicity,  we  drop  the  prime,  and  use  I(rj,  r,  x,/u)  for  /'(//,  r,  x,/u)  below. 


(Al) 


A.2.  The  WENO  algorithm:  approximations  to  the  spatial  derivatives 

The  spatial  derivative  terms  in  equation  (Al)  are  approximated  by  a  fifth-order  finite  differ¬ 
ence  WENO  scheme. 

We  first  give  the  WENO  reconstruction  procedure  in  approximating 

dl(r(\ruxj,nk )  1 


dx 


—(hj+i/2  ~  hj- 1/2) 


with  fixed  77  =  ij",  r  =  r,  and  /j  =  74  .  The  numerical  flux  hj+ 1/2  is  obtained  by  the  fifth-order  WENO 
approximation  in  an  upwind  fashion,  because  the  wind  direction  is  fixed  (negative).  Denote 

hj  =  I{rf‘ ,  ru  Xj,Hk),  j  =  -2,  -1,  •  •  • ,  Nx  +  3 

with  fixed  n,  i  and  k.  The  numerical  flux  from  the  WENO  procedure  is  obtained  by 

hj+ 1/2  =  kh^+1/2  +  +  a'3^)1/2,  (A2) 
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where  h('"\r  are  the  three  third-order  fluxes  on  three  different  stencils  given  by 


^/+ 1/2  =  + 

hj+\/2  ~  +  5^f+1  —  6^+2’ 

hj+ 1/2  =  ~  6^+2  +  3^7+3, 


and  the  nonlinear  weights  com  are  given  by 


(On 


0J„ 


D/= i  &/ 


oil 


7i 


(e  +  P,)2 


where  e  is  a  parameter  to  avoid  the  denominator  to  become  zero  and  is  taken  as  e  =  10  8.  The 
linear  weights  ji  are  given  by 


3  3  1 

7i  -  jo’  72  ~  5’  73  -  10’ 

and  the  smoothness  indicators  Pi  are  given  by 

13  ,1  , 

P\  -  —  (hj- 1  -  2/zy  +  hj+ 1 )  +  i  -  4/zy  +  3hj+i)  , 

p2  =  -  2/Z/+1  +  /27+2)2  +  -(/i,  -  hj+2)2, 

13  It 

p3  -  y^^O'+i  _  2/i7+2  +  hj+p)2  +  -(3/z/+i  -  4/t j+2  +  hj+p2. 


Similarly,  we  give  the  WENO  procedure  in  approximating 


d{l-p)l 

dii  ’ 


<9(1  -  n2)I(T]n,rhxk,Hj) 
d/j. 


^(^y+1/2  -  hj- 1/2) 


with  fixed  //  =  7/",  r  =  r,  and  .*  =  .vy .  The  numerical  flux  hj+\/2  is  also  obtained  by  the  fifth-order 
WENO  approximation  in  an  upwind  fashion,  however  the  wind  direction  here  is  positive,  opposite 
from  that  of  Denote 


hj  =  (1  -  n2)I(Tin,ri,xk,Hj), 


j  —  ~  3,  —2,  ••  •  ,Nfj  +  2 


with  fixed  n,  i  and  k.  The  numerical  flux  from  the  WENO  procedure  is  obtained  by 


hj+ 1/2  -  ^\h%r_  +  <y2^2,i/2  + 


(A3) 
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where  h('." '  are  the  three  third-order  fluxes  on  three  different  stencils  given  by 


0+1/2 


h{!\  -  --hj+2  +  —h j+\  +  -hj, 

O  O  3 


0+1/2 


h{2)  - 

2+i/“  3 


1,  5,  1, 

~hj+ 1  +  fhj  -  -hj- 1, 


/'(3)  -  11/,  .  _  1_U  .  +  -h  . 

y+ 1 /2  _  g  '0  g“/-l  +  3^2-2’ 


and  the  nonlinear  weights  cum  are  given  by 

_  6^/? 

— 


Z/=l  <22/ 


<22/  = 


7/ 


(6+ A)2 


where  £  is  taken  as  e  =  10  8.  The  linear  weights  y/  are  also  given  by 

3  3  1 

7l  “  jq.  72  -  y  73  -  jq, 

and  the  smoothness  indicators  /?/  are  given  by 

=  Y2  (^2+2  -  2/7 7+ 1  +  hj y  +  -(/ty+2  -  4/l  y+ 1  +  3hj)2, 
P2  =  I2^hj+l  ~  2hj  +  ^2-l)2  +  4(^+1  _  Vl)2’ 

A  =  J^(hj  -  2hj-\  +  hj-2)2  +  4(3 hj  -  4/z/_!  +  /7;-2)2. 


In  the  end,  we  approximate  the  r-derivative  in  equation  (Al),  following  the  reconstruction 
procedures  mentioned  above.  However,  we  need  to  check  the  wind  direction  at  the  r-boundary  of 
each  cell.  When  /j  >  0,  the  wind  direction  is  positive,  then  we  use  equation  (A3)  to  approximate 
the  numerical  flux,  while  when  /u  <  0,  we  use  equation  (A2). 


A.3.  High  order  numerical  integration 


The  integration  of  the  resonance  scattering  term  is  calculated  by  a  fifth  order  quadrature  (Shen 
et  al.2007) 


N 

right 

I  fip)diu  =  A//  2j  <2 >kf(Mk)  +  0(Aju5), 

deleft  k=l 
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where  Hk  -  Me  ft  +  (k  -  \)d/u  and  the  weights  are  defined  as, 


6463 

1457  741 

5537 

0)1  ~  5760’ 

0)1  ~  1920’  ^~  640’ 

“  5760’ 

_  5537 

741 

1457 

_  6463 

%""3  “  5760’ 

^~2  _  640’  “ 

1920’  ^ 

“  5760’ 

and  cok  =  1  otherwise. 


A.4.  Numerical  boundary  condition 

Following  Carrillo  et  al.  (2006),  at  /u  =  -1  and  /j  -  1 ,  we  take  the  boundary  conditions  as,  for 

d  >  o, 


I(ri,  r,x,-l-  n)  =  I(r],  r,  x,  -1  +  //), 

I(r/,  r,  x,  1  +  n)  =  I(tj,  r,  x,  1  -  /u), 

motivated  by  the  physical  meaning  of  //  as  the  cosine  of  the  angle  to  the  z-axis.  We  also  explicitly 
impose  h  i  =  hN/j+ 1  =  0  for  the  first  and  last  numerical  fluxes  in  order  to  enforce  conservation  of 
mass. 


A.5.  Time  evolution 

To  evolve  in  time,  we  use  the  third-order  TVD  Runge-Kutta  time  discretization  (Shu  &  Osher 
1988).  For  system  of  ODEs  ut  =  L(u),  the  third  order  Runge-Kutta  method  is 

«(1)  =  un  +  AtL(m",t"), 

u{2)  =  +  i(wa)  +  AtL(w(1),t''  +  Ar)), 

w',+1  =  +  |(M<2>  +  AtL(m(2),  t"  +  ^Ar)). 
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